Example use case: Binary black hole systems
In this notebook we will set up binary systems to find binary black hole systems. We will make use of the event based logs to find out whether a system has become a BBH and to store the summary of their evolution.
We run individual binary systems so we can play around with the input settings and see how that affects the formation of the BHBH system. There are other, more appropriate, ways to find populations of binary black hole systems. A good start for that is setting up a custom logging code and a parse function, and then run a grid of systems.
[1]:
import os
import json
from binarycpython.utils.functions import temp_dir, output_lines
from binarycpython.utils.run_system_wrapper import run_system
from binarycpython import Population
TMP_DIR = temp_dir("notebooks", "notebook_BHBH", clean_path=True)
EVENT_TYPE_INDEX = 3
VERBOSITY = 0
To run a system we can use the run_system function, which can parse function arguments for binary_c and will run a system. We will start with this but later we will use the evolve_single function from the Population object.
[2]:
# set filename
log_filename = os.path.join(TMP_DIR, "log_file.txt")
# Run via run_system
output = run_system(M_1=60,
M_2=40,
orbital_period=1,
BH_prescription='BH_BELCZYNSKI',
log_filename=log_filename,
wind_mass_loss='WIND_ALGORITHM_BINARY_C_2020',
clean_log='True', # removes log formatting
log_period_unit='day', # sets the orbital periods in the log output to days
api_log_filename_prefix=TMP_DIR)
# Readout the contents
with open(log_filename, 'r') as f:
print(f.read())
TIME M1/M☉ M2/M☉ st1 st2 SEP/R☉ PER ECC R1/ROL1 R2/ROL2 TYPE
RANDOM_SEED=13963 RANDOM_COUNT=0
0.0000 60.000 40.000 MS MS 19.536 1 0.00 1.503 1.384 "INITIAL "
0.0000 60.000 40.000 MS MS 19.536 1 0.00 1.503 1.384 Contact reached R/RL = 1.50329 1.38448, st 1 1, "Contact because R>RL"
0.0000 90.000 0.000 MS γ -1 -1 -1.00 0.000 0.000 "TYPE_CHNGE"
0.0000 90.000 0.000 MS γ -1 -1 -1.00 0.000 0.000 "BEG_RCHE 1>2"
0.0000 90.000 0.000 MS γ -1 -1 -1.00 0.000 0.000 "COALESCE"
0.0000 90.000 0.000 MS γ -1 -1 -1.00 0.000 0.000 "END_RCHE 1!>2"
6.6996 26.531 0.000 HG γ -1 -1 -1.00 0.000 0.000 "OFF_MS"
6.6996 26.531 0.000 HG γ -1 -1 -1.00 0.000 0.000 "TYPE_CHNGE"
6.7093 26.469 0.000 CHeB γ -1 -1 -1.00 0.000 0.000 "TYPE_CHNGE"
7.2750 9.883 0.000 HeMS γ -1 -1 -1.00 0.000 0.000 "TYPE_CHNGE"
7.5818 7.328 0.000 HeHG γ -1 -1 -1.00 0.000 0.000 "TYPE_CHNGE"
7.6220 2.909 0.000 BH γ -1 -1 -1.00 0.000 0.000 Randbuf=13963 - Mers(0)=0.0323703 - Mers(1)=0.973329 - Mers(2)=0.592129 - Mers(3)=0.0260156
7.6220 2.909 0.000 BH γ -1 -1 -1.00 0.000 0.000 SN kick Ib/c (SN type 13 13, pre-explosion M=6.92349 Mc"CO"=5.31476 stellar_type=8 Eexp=1 foe vexp=5004.45 km/s) -> kick 1(190) vk=255.754 vr=0 omega=0.163461 phi=0.185317 -> vn=255.754 ; final sep 0 ecc -1 (random count 0) - Runaway v=(252.345,40.9072,7.66879) |v|=255.754 : companion v=(0,0,0), |v|=0 ;
7.6220 2.909 0.000 BH γ -1 -1 -1.00 0.000 0.000 "TYPE_CHNGE"
7.6220 2.909 0.000 BH γ -1 -1 -1.00 0.000 0.000 "SN"
15000.0000 2.909 0.000 BH γ -1 -1 -1.00 0.000 0.000 "MAX_TIME"
From the output of this file we can see that the stellar types (K1, K2) are not entirely what we are looking for. One of them is γ, short for “massless remnant” i.e. no star any more, which means the stars must have merged to make a single star. To prevent this, we can star the system with a wider initial orbit, so that the stars do not interact (except by a little wind accretion).
[3]:
# run via run_system
output = run_system(M_1=60,
M_2=40,
orbital_period=2000, # days
BH_prescription='BH_BELCZYNSKI',
log_filename=log_filename,
wind_mass_loss='WIND_ALGORITHM_BINARY_C_2020',
clean_log='True', # removes log formatting
log_period_unit='day', # sets the orbital periods in the log output to days
api_log_filename_prefix=TMP_DIR)
# Readout contents
with open(log_filename, 'r') as f:
print(f.read())
TIME M1/M☉ M2/M☉ st1 st2 SEP/R☉ PER ECC R1/ROL1 R2/ROL2 TYPE
RANDOM_SEED=76943 RANDOM_COUNT=0
0.0000 60.000 40.000 MS MS 3101.2 2e+03 0.00 0.009 0.009 "INITIAL "
0.0000 60.000 40.000 MS MS 3101.2 2e+03 0.00 0.009 0.009 "BEG_SYMB"
4.3578 42.087 35.930 MS MS 3970.7 3.28e+03 0.00 0.023 0.014 "Start tidal lock 1"
4.3937 41.620 35.893 MS MS 3996.4 3.32e+03 0.00 0.023 0.014 "End tidal lock 1"
6.4207 27.619 27.607 MS MS 5602.5 6.53e+03 0.00 0.013 0.013 "Start tidal lock 2"
6.5391 27.162 27.156 MS MS 5695.9 6.75e+03 0.00 0.012 0.012 "End tidal lock 2"
6.6477 26.762 26.731 HG MS 5791.6 6.98e+03 0.00 0.010 0.011 "TYPE_CHNGE"
6.6550 26.729 26.728 HG MS 5794.4 6.99e+03 0.00 0.203 0.010 "Start tidal lock 2"
6.6552 26.727 26.728 HG MS 5794.6 6.99e+03 0.00 0.225 0.010 "q-inv"
6.6553 26.726 26.728 HG HG 5795.4 6.99e+03 0.00 0.231 0.010 "OFF_MS"
6.6553 26.726 26.728 HG HG 5795.4 6.99e+03 0.00 0.231 0.010 "TYPE_CHNGE"
6.6553 26.726 26.728 HG HG 5795.4 6.99e+03 0.00 0.231 0.010 "End tidal lock 2"
6.6556 26.722 26.728 HG HG 5796 6.99e+03 0.00 0.266 0.011 "Start tidal lock 2"
6.6559 26.717 26.727 HG HG 5796.5 6.99e+03 0.00 0.302 0.012 "End tidal lock 2"
6.6573 26.690 26.724 CHeB HG 5798.6 7e+03 0.00 0.546 0.023 "TYPE_CHNGE"
6.6649 26.505 26.682 CHeB CHeB 5812 7.04e+03 0.00 0.548 0.543 "TYPE_CHNGE"
7.2666 10.177 10.469 CHeB CHeB 8831.5 2.11e+04 0.00 0.000 0.458 "END_SYMB"
7.2673 10.171 10.449 HeMS CHeB 8840.7 2.12e+04 0.00 0.000 0.453 "TYPE_CHNGE"
7.2733 10.113 10.271 HeMS CHeB 8920.7 2.16e+04 0.00 0.000 0.263 "BEG_SYMB"
7.2804 10.031 10.162 HeMS HeMS 8982.4 2.19e+04 0.00 0.000 0.000 "TYPE_CHNGE"
7.4927 7.910 7.993 HeHG HeMS 11436 3.55e+04 0.00 0.000 0.000 "TYPE_CHNGE"
7.4988 7.854 7.943 HeHG HeHG 11513 3.6e+04 0.00 0.000 0.000 "TYPE_CHNGE"
7.5285 3.966 7.612 BH HeHG -23.473 -1 -1.00 0.000 0.000 Randbuf=76943 - Mers(0)=0.667522 - Mers(1)=0.671034 - Mers(2)=0.473751 - Mers(3)=0.362148 - Mers(4)=0.592519
7.5285 3.966 7.612 BH HeHG -23.473 -1 -1.00 0.000 0.000 SN kick Ib/c (SN type 13 13, pre-explosion M=7.49449 Mc"CO"=5.76446 stellar_type=8 Eexp=1 foe vexp=5338.47 km/s) -> kick 1(190) vk=320.854 vr=15.5127 omega=3.72291 phi=-0.279323 -> vn=307.408 ; final sep -23.4733 ecc -1 (random count 0) - Runaway v=(-260.125,-168.966,48.5359) |v|=313.959 : companion v=(-7.80471,-0.205277,0.0209149), |v|=7.80744 ;
7.5285 3.966 7.612 BH HeHG -23.473 -1 -1.00 0.000 0.000 "TYPE_CHNGE"
7.5285 3.966 7.612 BH HeHG -23.473 -1 -1.00 0.000 0.000 "DISRUPT "
7.5285 3.966 7.612 BH HeHG -23.473 -1 -1.00 0.000 0.000 "END_SYMB"
7.5285 3.966 7.612 BH HeHG -23.473 -1 -1.00 0.000 0.000 "SN"
7.5343 3.966 4.030 BH BH -23.473 -1 -1.00 0.000 0.000 Mers(5)=0.883215 - Mers(6)=0.0139536 - Mers(7)=0.79395 - Mers(8)=0.839891
7.5343 3.966 4.030 BH BH -23.473 -1 -1.00 0.000 0.000 SN kick Ib/c (SN type 13 13, pre-explosion M=7.52657 Mc"CO"=5.78976 stellar_type=8 Eexp=1 foe vexp=5362.7 km/s) -> kick 1(190) vk=433.264 vr=0 omega=5.27719 phi=0.62846 -> vn=433.264 ; final sep -23.4733 ecc -1 (random count 5) - Runaway v=(224.099,-296.256,-215.136) |v|=433.264 : companion v=(-6.54885e-51,-8.09379e-104,-2.38777e-104), |v|=6.54885e-51 ;
7.5343 3.966 4.030 BH BH -23.473 -1 -1.00 0.000 0.000 "TYPE_CHNGE"
7.5343 3.966 4.030 BH BH -23.473 -1 -1.00 0.000 0.000 "SN"
15000.0000 3.966 4.030 BH BH -23.473 -1 -1.00 0.000 0.000 "MAX_TIME"
It would be great if we can automatically detect if a system becomes a BH-BH system. We can do this with the event based logging. One of the events that is available is the DCO_formation, which stands for ‘double compact-object formation’. If this event is present in the output then we know a compact object of some sort was formed. If we then extract the correct parameters we can find out if there is a BHBH system.
Lets set up a population object and configure it. We then want to use the evolve_single function to run the systems.
[4]:
BHBH_pop = Population(verbosity=VERBOSITY, tmp_dir=TMP_DIR)
# binary-c options related to event-based logging
BHBH_pop.set(
event_based_logging_RLOF=1,
event_based_logging_DCO=1,
)
# Configure system
BHBH_pop.set(
M_1=60,
M_2=40,
orbital_period=2000, # days
BH_prescription='BH_BELCZYNSKI',
log_filename=log_filename,
wind_mass_loss='WIND_ALGORITHM_BINARY_C_2020',
clean_log='True', # removes log formatting
log_period_unit='day', # sets the orbital periods in the log output to days
)
[5]:
BHBH_pop.evolve_single()
[5]:
'SINGLE_STAR_LIFETIME 60 7.52846\nEVENT C02240E9-6E49-4938-9E9D-9C1DEE09400B 1 0 DCO_formation 60 40 2000 3101.19 0 7.534313670863e+00 0.02 58237 14 14 3.96639 4.03029 -4.47113 -1 106.889 -4.47113 -1 106.889 7.53431e+06 4.69096e+08 4.7663e+08 0 0 0\n'
The formatting of this output is not ideal, so lets write a functions to extract the events and parse the info. We can use the event header dictionaries to parse the events properly.
[6]:
event_based_logging_parameter_list_dict = BHBH_pop.population_options['event_based_logging_parameter_list_dict']
def recast_values(event_dict):
"""
Function to recast the values from strings to number values
"""
for key in event_dict.keys():
if key in ['uuid', 'event_type']:
continue
try:
if '.' in event_dict[key]:
event_dict[key] = float(event_dict[key])
else:
event_dict[key] = int(event_dict[key])
except:
event_dict[key] = float(event_dict[key])
return event_dict
def parse_output_evolution(evolution_output, parsing_dict):
"""
Function to parse the output of the evolution of the system and create a dictionary containing the
"""
events_list = []
# Loop over output
for line in output_lines(evolution_output):
if line.startswith("EVENT"):
# Parse output and create dictionary
parameter_values = line.split()[1:]
event_type = parameter_values[EVENT_TYPE_INDEX]
parameter_names = parsing_dict[event_type]
event_dict = {parameter_name: parameter_value for (parameter_name, parameter_value) in zip(parameter_names, parameter_values)}
# recast values
event_dict = recast_values(event_dict=event_dict)
#
events_list.append(event_dict)
return events_list
def find_BHBH_formation_event(events_list):
"""
Function to extract a BHBH formation event from the events list
"""
found_BHBH_formation_event = False
BHBH_formation_event = {}
for event in events_list:
if (event['event_type'] == "DCO_formation") and (event['DCO_stellar_type_1'] == 14 and event['DCO_stellar_type_2']==14):
found_BHBH_formation_event = True
BHBH_formation_event = event
return found_BHBH_formation_event, BHBH_formation_event
# Evolve
evolution_output = BHBH_pop.evolve_single()
# parse output
events_list = parse_output_evolution(evolution_output=evolution_output, parsing_dict=event_based_logging_parameter_list_dict)
# Extract BHBH formation event
found_BHBH_formation_event, BHBH_formation_event = find_BHBH_formation_event(events_list=events_list)
# Print info
if not found_BHBH_formation_event:
print("No BHBH formation event found")
else:
if BHBH_formation_event['DCO_separation'] < 0:
print("BHBH system is unbound")
else:
print("BHBH system is bound")
BHBH system is unbound
Now that we have a function to detect a BHBH system, we want to set up a search for these systems.
[7]:
def search_for_BHBH(BHBH_pop, maxcount, verbosity=0):
"""
Function to search for BHBH formation events
"""
#
found_bound_BHBH_system = False
count = 0
#
while found_bound_BHBH_system == False and count < maxcount:
count = count + 1
if verbosity: print("system {} / {}".format(count,maxcount))
# Evolve
evolution_output = BHBH_pop.evolve_single()
# parse output
events_list = parse_output_evolution(evolution_output=evolution_output, parsing_dict=event_based_logging_parameter_list_dict)
# Extract BHBH formation event
found_BHBH_formation_event, BHBH_formation_event = find_BHBH_formation_event(events_list=events_list)
# Check if we have a bound BHBH system
if found_BHBH_formation_event:
if BHBH_formation_event['DCO_separation'] > 0.0:
print('Found bound BHBH system')
found_bound_BHBH_system = True
return BHBH_formation_event, events_list
# If we did not find any BHBH system we return an empty dict
if found_bound_BHBH_system == False:
print("No BHBH systems found")
return {}, []
[8]:
# Configure population
BHBH_pop.set(
M_1=60,
M_2=40,
BH_prescription='BH_BELCZYNSKI',
orbital_period=2000, # days
wind_mass_loss='WIND_ALGORITHM_BINARY_C_2020',
clean_log='True', # removes log formatting
log_period_unit='day', # sets the orbital periods in the log output to days
)
number_of_systems = 100
search_result = search_for_BHBH(
BHBH_pop=BHBH_pop,
maxcount=number_of_systems,
verbosity=VERBOSITY,
)
print(search_result)
No BHBH systems found
({}, [])
How can we help the system become a BHBH merger? We can try turning off SN kicks.
[9]:
# Configure population
BHBH_pop.set(
sn_kick_dispersion_II=0,
sn_kick_dispersion_IBC=0,
sn_kick_dispersion_GRB_COLLAPSAR=0,
)
number_of_systems = 100
search_result = search_for_BHBH(
BHBH_pop=BHBH_pop,
maxcount=number_of_systems,
verbosity=VERBOSITY,
)
print(search_result)
Found bound BHBH system
({'uuid': '74225838-9B45-4A27-AACA-92AF2508D184', 'probability': 1, 'event_number': 0, 'event_type': 'DCO_formation', 'zams_mass_1': 60, 'zams_mass_2': 40, 'zams_orbital_period': 2000, 'zams_separation': 3101.19, 'zams_eccentricity': 0, 'time': 7.534313670863, 'metallicity': 0.02, 'random_seed': 17423, 'DCO_stellar_type_1': 14, 'DCO_stellar_type_2': 14, 'DCO_mass_1': 3.96639, 'DCO_mass_2': 4.03029, 'DCO_separation': 90146.7, 'DCO_eccentricity': 0.865727, 'DCO_period': 3035.13, 'DCO_previous_separation': 17263.9, 'DCO_previous_eccentricity': 0.304709, 'DCO_previous_period': 211.84, 'DCO_formation_time_in_years': 7534310.0, 'DCO_inspiral_time_in_years': 5.46768e+23, 'DCO_merger_time_in_years': 5.46768e+23, 'DCO_total_rlof_episodes': 0, 'DCO_stable_rlof_episodes': 0, 'DCO_unstable_rlof_episodes': 0}, [{'uuid': '74225838-9B45-4A27-AACA-92AF2508D184', 'probability': 1, 'event_number': 0, 'event_type': 'DCO_formation', 'zams_mass_1': 60, 'zams_mass_2': 40, 'zams_orbital_period': 2000, 'zams_separation': 3101.19, 'zams_eccentricity': 0, 'time': 7.534313670863, 'metallicity': 0.02, 'random_seed': 17423, 'DCO_stellar_type_1': 14, 'DCO_stellar_type_2': 14, 'DCO_mass_1': 3.96639, 'DCO_mass_2': 4.03029, 'DCO_separation': 90146.7, 'DCO_eccentricity': 0.865727, 'DCO_period': 3035.13, 'DCO_previous_separation': 17263.9, 'DCO_previous_eccentricity': 0.304709, 'DCO_previous_period': 211.84, 'DCO_formation_time_in_years': 7534310.0, 'DCO_inspiral_time_in_years': 5.46768e+23, 'DCO_merger_time_in_years': 5.46768e+23, 'DCO_total_rlof_episodes': 0, 'DCO_stable_rlof_episodes': 0, 'DCO_unstable_rlof_episodes': 0}])
You should now have found a BHBH system but usually they have very wide orbits. This is caused by mass loss before, and during, the supernova.
We can reduce the former by setting the wind mass loss to zero. This can be done unphysically by setting wind_mass_loss=0, or you can reduce the massive-star winds by setting the metallicity to be small, e.g. \(10^{-3}\).
[10]:
BHBH_pop.set(
metallicity=0.001,
)
number_of_systems = 100
BHBH_formation_event, events_list = search_for_BHBH(
BHBH_pop=BHBH_pop,
maxcount=number_of_systems,
verbosity=VERBOSITY,
)
for event in events_list:
print(event['event_type'], "\n", event, "\n\n")
Found bound BHBH system
RLOF
{'uuid': '7B34580B-BA86-4364-B219-7BE72F515970', 'probability': 1, 'event_number': 0, 'event_type': 'RLOF', 'zams_mass_1': 60, 'zams_mass_2': 40, 'zams_orbital_period': 2000, 'zams_separation': 3101.19, 'zams_eccentricity': 0, 'time': 4.058557979452, 'metallicity': 0.001, 'random_seed': 36879, 'RLOF_initial_mass_accretor': 39.7575, 'RLOF_initial_mass_donor': 54.5101, 'RLOF_initial_radius_accretor': 11.3814, 'RLOF_initial_radius_donor': 1319.94, 'RLOF_initial_separation': 3246.27, 'RLOF_initial_orbital_period': 6.04091, 'RLOF_initial_stellar_type_accretor': 1, 'RLOF_initial_stellar_type_donor': 4, 'RLOF_initial_orbital_angular_momentum': 251989000.0, 'RLOF_initial_stability': 0, 'RLOF_initial_starnum_accretor': 1, 'RLOF_initial_starnum_donor': 0, 'RLOF_initial_time': 3890959.062818, 'RLOF_initial_disk': 1, 'RLOF_final_mass_accretor': 50.4535, 'RLOF_final_mass_donor': 24.9596, 'RLOF_final_radius_accretor': 10.547, 'RLOF_final_radius_donor': 1012.85, 'RLOF_final_separation': 3211.47, 'RLOF_final_orbital_period': 6.64565, 'RLOF_final_stellar_type_accretor': 1, 'RLOF_final_stellar_type_donor': 4, 'RLOF_final_orbital_angular_momentum': 162835000.0, 'RLOF_final_stability': 0, 'RLOF_final_starnum_accretor': 1, 'RLOF_final_starnum_donor': 0, 'RLOF_final_time': 4058557.979452, 'RLOF_final_disk': 0, 'RLOF_total_mass_lost': 0, 'RLOF_total_mass_accreted': 4.41063, 'RLOF_total_mass_transferred': 4.41063, 'RLOF_total_mass_lost_from_accretor': 0, 'RLOF_total_mass_lost_from_common_envelope': 0, 'RLOF_total_time_spent_masstransfer': 167590, 'RLOF_episode_number': 1}
RLOF
{'uuid': '7B34580B-BA86-4364-B219-7BE72F515970', 'probability': 1, 'event_number': 1, 'event_type': 'RLOF', 'zams_mass_1': 60, 'zams_mass_2': 40, 'zams_orbital_period': 2000, 'zams_separation': 3101.19, 'zams_eccentricity': 0, 'time': 6.10177351634, 'metallicity': 0.001, 'random_seed': 36879, 'RLOF_initial_mass_accretor': 22.2425, 'RLOF_initial_mass_donor': 39.1412, 'RLOF_initial_radius_accretor': 9.43081e-05, 'RLOF_initial_radius_donor': 1669.7, 'RLOF_initial_separation': 3894.32, 'RLOF_initial_orbital_period': 9.83617, 'RLOF_initial_stellar_type_accretor': 14, 'RLOF_initial_stellar_type_donor': 4, 'RLOF_initial_orbital_angular_momentum': 137399000.0, 'RLOF_initial_stability': 0, 'RLOF_initial_starnum_accretor': 0, 'RLOF_initial_starnum_donor': 1, 'RLOF_initial_time': 5978661.25158, 'RLOF_initial_disk': 1, 'RLOF_final_mass_accretor': 22.2562, 'RLOF_final_mass_donor': 20.582, 'RLOF_final_radius_accretor': 9.43663e-05, 'RLOF_final_radius_donor': 1683.78, 'RLOF_final_separation': 4541.55, 'RLOF_final_orbital_period': 14.8282, 'RLOF_final_stellar_type_accretor': 14, 'RLOF_final_stellar_type_donor': 4, 'RLOF_final_orbital_angular_momentum': 93459900.0, 'RLOF_final_stability': 0, 'RLOF_final_starnum_accretor': 0, 'RLOF_final_starnum_donor': 1, 'RLOF_final_time': 6101773.51634, 'RLOF_final_disk': 0, 'RLOF_total_mass_lost': -2.89615, 'RLOF_total_mass_accreted': 0.0923439, 'RLOF_total_mass_transferred': 0.0923439, 'RLOF_total_mass_lost_from_accretor': -2.89615, 'RLOF_total_mass_lost_from_common_envelope': 0, 'RLOF_total_time_spent_masstransfer': 123101, 'RLOF_episode_number': 2}
DCO_formation
{'uuid': '7B34580B-BA86-4364-B219-7BE72F515970', 'probability': 1, 'event_number': 2, 'event_type': 'DCO_formation', 'zams_mass_1': 60, 'zams_mass_2': 40, 'zams_orbital_period': 2000, 'zams_separation': 3101.19, 'zams_eccentricity': 0, 'time': 6.337299289051, 'metallicity': 0.001, 'random_seed': 36879, 'DCO_stellar_type_1': 14, 'DCO_stellar_type_2': 14, 'DCO_mass_1': 22.2564, 'DCO_mass_2': 19.9264, 'DCO_separation': 4597.19, 'DCO_eccentricity': 0, 'DCO_period': 15.2187, 'DCO_previous_separation': 4597.17, 'DCO_previous_eccentricity': 0, 'DCO_previous_period': 15.2186, 'DCO_formation_time_in_years': 6337300.0, 'DCO_inspiral_time_in_years': 3.58246e+18, 'DCO_merger_time_in_years': 3.58246e+18, 'DCO_total_rlof_episodes': 2, 'DCO_stable_rlof_episodes': 2, 'DCO_unstable_rlof_episodes': 0}
Oh dear! While the mass loss is reduced, there is now Roche-lobe overflow (RLOF) in the system. Two stable RLOF episodes even!
The first two events in the events_list are RLOF event types. The first RLOF event is RLOF from the primary ('RLOF_initial_starnum_donor': 0) to the secondary 'RLOF_initial_starnum_accretor': 1).
This increases the secondary’s mass ('RLOF_initial_mass_accretor': 39.7575 and 'RLOF_final_mass_accretor': 50.4535,)
We can try making the system wider to prevent this, but we want a shorter period! Instead, let’s shorten the initial period to force common-envelope evolution. This will shrink the system. We will also turn off the mass loss to give us the best change of acquiring the system we want, and turn the SN kicks back on (they have less effect in closer, more gravitationally-bound systems).
[11]:
BHBH_pop.set(
metallicity=0.0001,
orbital_period=2, # days
wind_mass_loss='WIND_ALGORITHM_NONE',
alpha_ce = 1,
lambda_ce = 0.5,
sn_kick_dispersion_II=190,
sn_kick_dispersion_IBC=190,
sn_kick_dispersion_GRB_COLLAPSAR=190
)
number_of_systems = 10
BHBH_formation_event, events_list = search_for_BHBH(
BHBH_pop=BHBH_pop,
maxcount=number_of_systems,
verbosity=VERBOSITY,
)
for event in events_list:
print(event['event_type'], "\n", event, "\n\n")
Found bound BHBH system
RLOF
{'uuid': 'A4E35A70-0A50-485E-954D-82264F485D96', 'probability': 1, 'event_number': 0, 'event_type': 'RLOF', 'zams_mass_1': 60, 'zams_mass_2': 40, 'zams_orbital_period': 2, 'zams_separation': 31.0119, 'zams_eccentricity': 0, 'time': 3.918234836492, 'metallicity': 0.0001, 'random_seed': 4892, 'RLOF_initial_mass_accretor': 40, 'RLOF_initial_mass_donor': 60, 'RLOF_initial_radius_accretor': 7.88211, 'RLOF_initial_radius_donor': 12.5742, 'RLOF_initial_separation': 30.3393, 'RLOF_initial_orbital_period': 0.00529926, 'RLOF_initial_stellar_type_accretor': 1, 'RLOF_initial_stellar_type_donor': 1, 'RLOF_initial_orbital_angular_momentum': 26193100.0, 'RLOF_initial_stability': 0, 'RLOF_initial_starnum_accretor': 1, 'RLOF_initial_starnum_donor': 0, 'RLOF_initial_time': 2895212.520129, 'RLOF_initial_disk': 0, 'RLOF_final_mass_accretor': 77.16, 'RLOF_final_mass_donor': 22.84, 'RLOF_final_radius_accretor': 9.6474, 'RLOF_final_radius_donor': 12.9501, 'RLOF_final_separation': 52.9579, 'RLOF_final_orbital_period': 0.0122209, 'RLOF_final_stellar_type_accretor': 1, 'RLOF_final_stellar_type_donor': 4, 'RLOF_final_orbital_angular_momentum': 25411300.0, 'RLOF_final_stability': 0, 'RLOF_final_starnum_accretor': 1, 'RLOF_final_starnum_donor': 0, 'RLOF_final_time': 3918234.836492, 'RLOF_final_disk': 0, 'RLOF_total_mass_lost': 0, 'RLOF_total_mass_accreted': 37.16, 'RLOF_total_mass_transferred': 37.16, 'RLOF_total_mass_lost_from_accretor': 0, 'RLOF_total_mass_lost_from_common_envelope': 0, 'RLOF_total_time_spent_masstransfer': 1021860.0, 'RLOF_episode_number': 1}
RLOF
{'uuid': 'A4E35A70-0A50-485E-954D-82264F485D96', 'probability': 1, 'event_number': 1, 'event_type': 'RLOF', 'zams_mass_1': 60, 'zams_mass_2': 40, 'zams_orbital_period': 2, 'zams_separation': 31.0119, 'zams_eccentricity': 0, 'time': 6.047398220502, 'metallicity': 0.0001, 'random_seed': 4892, 'RLOF_initial_mass_accretor': 22.8955, 'RLOF_initial_mass_donor': 75.2163, 'RLOF_initial_radius_accretor': 9.70767e-05, 'RLOF_initial_radius_donor': 203.37, 'RLOF_initial_separation': 35.0383, 'RLOF_initial_orbital_period': 0.00663959, 'RLOF_initial_stellar_type_accretor': 14, 'RLOF_initial_stellar_type_donor': 4, 'RLOF_initial_orbital_angular_momentum': 20392800.0, 'RLOF_initial_stability': 3, 'RLOF_initial_starnum_accretor': 0, 'RLOF_initial_starnum_donor': 1, 'RLOF_initial_time': 6044477.100382, 'RLOF_initial_disk': 1, 'RLOF_final_mass_accretor': 22.8955, 'RLOF_final_mass_donor': 31.9069, 'RLOF_final_radius_accretor': 9.70767e-05, 'RLOF_final_radius_donor': 1.90098, 'RLOF_final_separation': 8.56607, 'RLOF_final_orbital_period': 0.00107394, 'RLOF_final_stellar_type_accretor': 14, 'RLOF_final_stellar_type_donor': 7, 'RLOF_final_orbital_angular_momentum': 5722660.0, 'RLOF_final_stability': 3, 'RLOF_final_starnum_accretor': 0, 'RLOF_final_starnum_donor': 1, 'RLOF_final_time': 6047398.220502, 'RLOF_final_disk': 0, 'RLOF_total_mass_lost': 41.4211, 'RLOF_total_mass_accreted': 1.94368, 'RLOF_total_mass_transferred': 45.2531, 'RLOF_total_mass_lost_from_accretor': -1.88822, 'RLOF_total_mass_lost_from_common_envelope': 43.3094, 'RLOF_total_time_spent_masstransfer': 637470, 'RLOF_episode_number': 3}
DCO_formation
{'uuid': 'A4E35A70-0A50-485E-954D-82264F485D96', 'probability': 1, 'event_number': 2, 'event_type': 'DCO_formation', 'zams_mass_1': 60, 'zams_mass_2': 40, 'zams_orbital_period': 2, 'zams_separation': 31.0119, 'zams_eccentricity': 0, 'time': 6.375786951085, 'metallicity': 0.0001, 'random_seed': 4892, 'DCO_stellar_type_1': 14, 'DCO_stellar_type_2': 14, 'DCO_mass_1': 22.8955, 'DCO_mass_2': 31.9069, 'DCO_separation': 16.5045, 'DCO_eccentricity': 0.494521, 'DCO_period': 0.00287219, 'DCO_previous_separation': 8.34723, 'DCO_previous_eccentricity': 1.41494e-06, 'DCO_previous_period': 0.00103305, 'DCO_formation_time_in_years': 6375790.0, 'DCO_inspiral_time_in_years': 95795100.0, 'DCO_merger_time_in_years': 102171000.0, 'DCO_total_rlof_episodes': 3, 'DCO_stable_rlof_episodes': 2, 'DCO_unstable_rlof_episodes': 1}
Note that this system has a far shorter period. Let’s run a number of these to see what we find.
Lets now try to find a system which at the formation of the final black hole has a separation of 30 \(R_{\odot}\) or less
[12]:
found = False
while found == False:
number_of_systems = 10
BHBH_formation_event, events_list = search_for_BHBH(
BHBH_pop=BHBH_pop,
maxcount=number_of_systems,
verbosity=VERBOSITY,
)
if events_list and BHBH_formation_event['DCO_separation'] < 30:
found = True
print(json.dumps(BHBH_formation_event, indent=4))
Found bound BHBH system
{
"uuid": "145CDED3-2593-4483-BEB8-D1A1BE5D9018",
"probability": 1,
"event_number": 2,
"event_type": "DCO_formation",
"zams_mass_1": 60,
"zams_mass_2": 40,
"zams_orbital_period": 2,
"zams_separation": 31.0119,
"zams_eccentricity": 0,
"time": 6.370900869438,
"metallicity": 0.0001,
"random_seed": 22112,
"DCO_stellar_type_1": 14,
"DCO_stellar_type_2": 14,
"DCO_mass_1": 22.8401,
"DCO_mass_2": 32.5402,
"DCO_separation": 7.0094,
"DCO_eccentricity": 0.40487,
"DCO_period": 0.000790771,
"DCO_previous_separation": 9.84471,
"DCO_previous_eccentricity": 6.5446e-06,
"DCO_previous_period": 0.00131624,
"DCO_formation_time_in_years": 6370900.0,
"DCO_inspiral_time_in_years": 4414940.0,
"DCO_merger_time_in_years": 10785800.0,
"DCO_total_rlof_episodes": 3,
"DCO_stable_rlof_episodes": 2,
"DCO_unstable_rlof_episodes": 1
}
Playing around with the input for an individual system gives us some understanding of which parameters affect the formation of black hole systems, but, as mentioned in the introduction, this is not the most appropriate way to find many such systems. We want to make use of the Population utilities, and set up custom logging and a parse function.
Things to try next:
Set up a population object and set up sampling in M1, q and orbital period (e.g. 10 x 10 x 10)
Evolve this population at various metallicities
Capture all the BHBH DCO_formation events
Do you see a trend?